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Tine Fermi Gamma-ray Space Telescope has revealed a diffuse y-ray 
ackground at energies from 0.1 GeV to 1 TeV, which can be sepa- 
ted into Galactic emission and an isotropic, extragalactic compo- 
Ghent 1, Previous efforts to understand the latter have been hampered 
Hy the lack of physical models capable of predicting the y-ray emis- 
sion produced by the many candidate sources, primarily active galac- 
ic nuclei 7° and star-forming galaxies f", leaving their contributions 
poorly constrained. Here we present a calculation of the contribution 
—of star-forming galaxies to the y-ray background that does not rely on 
empirical scalings, and is instead based on a physical model for the 
O@-ray emission produced when cosmic rays accelerated in supernova 
COvyemnants interact with the interstellar medium “!. After validating the 
Wmodel against local observations, we apply it to the observed cosmo- 
(“iogical star-forming galaxy population and recover an excellent match 
to both the total intensity and the spectral slope of the y-ray back- 
ground, demonstrating that star-forming galaxies alone can explain 
Che full diffuse, isotropic y-ray background. 


Many candidate sources have been proposed for the origin of the dif- 
Chie isotropic y-ray background. These include active galactic nuclei 
AGN) (particularly blazars 7), millisecond pulsars *, star-forming galax- 
bes (SFGs) 6-10 and dark matter annihilation '*. Previous estimates of their 
ontributions have relied on a highly-uncertain process of empirically scal- 
caing the emission from a small sample of local, resolved sources by their 
estimated cosmological abundances, whereas our approach in this paper 
is instead to calculate the emission from SFGs directly using a physical 
model. The cosmic rays (CRs) responsible for y-ray emission in SFGs 
(including the Milky Way) are produced by diffusive acceleration at su- 
pernova remnant shocks "°. This process transfers ~10% of the supernova 
mechanical energy to relativistic ions, yielding on average ~ 10°° erg in 
CR ions per supernova '*'>, with another ~2% (~ 2 x 10*° erg) deposited 
in CR electrons '®. The resulting CRs follow a power law distribution in 
particle momentum p of the form dn/dp x p~? 17.18. observations of indi- 
vidual supernova remnants, analytical models, and numerical simulations 
all indicate that the index q is in the range q ~ 2.0 — 2.6, with a mean 
value of q ~ 2.2 — 2.3 °°. Some of the CR ions collide inelastically 
with interstellar medium (ISM) nuclei, producing roughly equal numbers 
of x°, n and m~ mesons that rapidly decay via the channels 7° —> 2y, 
nm > u` +0, andat + u* +v,. The decay of 7° particles is respon- 


sible for most of the observed Galactic y-ray foreground, which displays a 
characteristic spectrum that rises sharply from ~ 0.1 — 1 GeV as a result 
of the 135 MeV rest mass of the 7° particle. 


As this discussion suggests, a SFG’s diffuse y-ray emission depends 
primarily on three factors: its total star formation rate (which determines 
its supernova rate and thus the rate at which CRs are injected), the distri- 
bution of y-ray energies produced when individual CRs collide with ISM 
nuclei (which depends on the parent CR energy E), and the fraction of CRs 
(again as a function of E) that undergo inelastic collisions before escaping 
the galaxy. The first two of these are relatively well-understood, but the 
third factor, known as the calorimetry fraction fcai(), is much less cer- 
tain. It depends on the properties of the galaxy, and the lack of a model 
for this dependence has previously precluded direct calculation of SFG y- 
ray emission. However, Ref. u recently introduced a model for fcai(F), 
based on rates of CR diffusion determined by the balance between the CR 
streaming instability and ion-neutral damping. In the Methods we describe 
a new technique to use this model to compute f-a1(/), and thus the total 
y-ray emission produced by CR ions in SFGs. 


We supplement the y-ray production rate from CR ions by adding the 
contribution from both primary CR electrons, directly injected by super- 
nova remnants, and secondary CR leptons (electrons and positrons), pro- 
duced in the 7* decay chain; these become important at energies < 1 
GeV. Our model for these particles includes energy losses due to ionisa- 
tion, synchrotron emission, bremsstrahlung, and inverse Compton scatter- 
ing. The model also includes the attenuation of y-rays produced by both 
CR ions and leptons due to pair production in collisions with far-infrared 
photons inside the source galaxy and extragalactic background light pho- 
tons outside the galaxy, which become important at energies > 100 GeV. 
The radiation that is absorbed by the host galaxy and extragalactic pho- 
ton fields is reprocessed to lower energies in a pair-production cascade, 
whereby the initial high energy pair inverse Compton scatters lower en- 
ergy photons up to y-ray energies and these, in turn, produce further pairs, 
and so on. Details of our calculation of all these processes are provided in 
the Methods. 


We now have a model that predicts the y-ray emission of a SFG. The 
next step in our analysis is to apply this model to a galaxy survey that 
samples the SFG population out to the epoch of peak cosmological star 
formation at z ~ 2. For this purpose, we make use of the Cosmic Assem- 
bly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) 7!” in 
the GOODS-S field. We apply our model to the CANDELS sample as 
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Figure 1: The y-ray spectra of nearby SFGs Predicted (lines) and ob- 
served (points) spectra for a selection of nearby SFGs detected in y-rays. 
The observations shown are taken from a combination of Fermi LAT ”, 
HESS ’, and VERITAS ™ where the horizontal bars show the energy bin 
and the vertical bars the 1 o uncertainty limit; in Panel a we show the 
local starburst galaxies Arp 220 and NGC 253, and Panel b we show the 
local quiescent galaxies M31 and NGC 4945. The solid lines show model 
predictions using only stellar data of the type we have available for the 
CANDELS sample, while the dotted lines shows results predicted if we 
supplement this with observed gas data. We list the full set of observed 
quantities used in computing these models in Extended Data Table 1. 


described in the Methods. To verify that our approach predicts reason- 
ably accurate y-ray spectra, we apply it to four local, resolved galaxies 
with measured y-ray emission ’7**, chosen to span a wide range of gas 
and star formation surface densities: Arp 220, NGC 253, M31, and NGC 
4945. The input data we use for these calculations are summarised in Ex- 
tended Data Table 1, and we show the results of the computation in Figure 
1, where the solid lines show the spectra derived using only stellar data 
(as we have for CANDELS) and, for comparison, the dotted lines show 
the results we obtain if we add directly-measured gas properties (available 
for these local galaxies). We see that the fits are slightly improved if we 
make direct use of gas data but, even for the stellar data only, our model 
reproduces the observed y-ray spectra to better than a factor of 2 for all 
galaxies at energies > 1 GeV, and within a factor of ~ 1.5 for the two 
more rapidly star-forming galaxies, which, as we show below, are more 
akin to the population that dominates the y-ray background. 

Having verified that we can obtain accurate predictions of y-ray spec- 
tra from stellar data alone, we carry out two additional validation steps. 
First, we examine the correlation between galaxies’ far-infrared and to- 
tal y-ray luminosities, computed as described in the Methods. In Fig- 
ure 2 we show the resulting distribution of galaxies in the Lr - Ly 
plane, along with a power law fit to the data (blue line), L,/ergs~+ = 
1028-2641.55 (Lrir/Lo)" 14+016, Our model prediction shows good 
agreement with the observed relation 7°, and we note that both the model 
and the observed correlation differ noticeably from the calorimetric limit 
obtained by simply setting fcai (E) = 1 (red line in the Figure). Thus 
the agreement is non-trivial, and suggests that our model is correctly pre- 
dicting the variation in galaxies’ calorimetry fractions as a function of star 
formation rate. 


Our second validation test is to compare our model with counts of re- 
solved SFGs observed by Fermi LAT. Details of how we perform the com- 
parison are given in the Methods. We show the results in Figure 3, which 
demonstrates that our model predicts SFG source counts consistent with 
observations, with the exception that we do not predict sources as bright 
as the Milky Way’s two satellite galaxies, the Large and Small Magellanic 
Clouds. This is not surprising, since our comparison includes only field 
galaxies. 
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Figure 2: The FIR-y correlation The correlation between far-infrared 
(8 — 1000 um) and y-ray (0.1 — 100 GeV) luminosity for the CAN- 
DELS sample, derived using our model. Points show individual CAN- 
DELS galaxies, colour-coded by redshift z. The blue line is a power law 
fit to the CANDELS sample with the shaded band containing 90% of data 
around the model fit. For comparison, the solid green line shows the em- 
pirical relation measured for 14 nearby, resolved SFGs * with 2 ø uncer- 
tainty in the shaded band. The red line is the calorimetric limit obtained 
by taking fcai = 1 at all energies in{Equation 1] as obtained by Ref. 7. 


Our final step is to compute the contribution of SFGs to the diffuse, 
isotropic y-ray background (see Methods for details). We present the re- 
sults of this calculation in Figure 4, and provide a detailed analysis of the 
model uncertainties in the Supplementary Information and Extended Data 
Figure 1. Figure 4 shows that the expected contribution of SFGs to the dif- 
fuse isotropic y-ray background fully reproduces both the intensity and the 
spectral shape of the observations from ~ 0.2 GeV to ~ 1 TeV. We empha- 
sise that we obtain this agreement from our model with no free parameters: 
our only inputs are the CR injection spectral index (q = 2.2), the energy 
per supernova (10°* erg), and the fraction of supernova energy that goes 
into primary ions and electrons (10% and 2%, respectively) — all quantities 
that are directly measured in the local Universe — and the distribution of 
SFGs sampled by CANDELS. The key to the success of the model is the 
galaxy-by-galaxy calculation of the energy-dependent calorimetry fraction 
fcai(E), which we demonstrate by also plotting the result (dotted line) we 
would obtain simply by setting fca1 = 1 for all galaxies at all energies. 
This clearly both overestimates the intensity and yields a spectral slope 
that is flatter than observed. 

We show the relative contributions to the background made by galaxies 
with differing star formation rates and redshifts in Extended Data Figure 
2. The Figure shows that the background at lower energies is dominated 
by galaxies from just after cosmic noon (z ~ 1 — 2), while at higher 
energies, where attenuation by extragalactic background light has a larger 
effect, the dominant contribution shifts towards lower redshifts, so that at 1 
TeV the background is dominated by z ~ 0.1 sources. At all energies, the 
dominant contribution comes from galaxies at the upper end of the star- 
forming main sequence, which have high but not extreme star formation 
rates for their redshift. 

It is important to put our finding that SFGs dominate the diffuse, 
isotropic y-ray background in the context of recent work, where a number 
of authors have argued that blazars and other AGN sources contribute sub- 
stantially or even dominate the background. We provide a more detailed 
discussion in the Supplementary Information, but here note that we find 
that, while blazars dominate the resolved component of the extragalactic 
y-ray background, as shown in Figure 3, SFGs dominate the unresolved 
component. This finding is consistent with statistical analyses of angular 
fluctuations in the isotropic background and cross-correlations between it 
and galaxies and quasars, which strongly disfavour blazars as a dominant 
contributor **°?°, Indeed, a straightforward extrapolation of the number 
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Figure 3: The y-ray source count distribution Red points show the compensated distribution of SFG luminosities > S?(dN/dS) on the sky as a 
function of photon flux S integrated from 1 - 100 GeV as seen by Fermi LAT; error bars show 90% confidence intervals or upper limits. The three 
brightest non-empty bins each contain only a single SFG, which we have labelled. Green points show model-predicted source counts for observed 
CANDELS galaxies at z > 0.1 with the 90% confidence limit, and blue points show Monte Carlo realisations of the z < 0.1 SFG population, with the 
light and dark shaded bands indicating 68% and 90% confidence intervals. Black squares show Fermi-detected blazars, and the orange band shows the 
blazar distribution model of Ref. 7’ within the 1 ø range. Finally, the red vertical band indicates the flux range over which Fermi observations become 
incomplete; the left edge of this band is the 4FGL threshold for 98% detection efficiency for sources with spectral index 2.3 7’. 


also suggests that blazars do not dominate the unresolved background. Our 
finding that SFGs alone are able to reproduce the full background is also 
consistent with the conclusions of Refs. °? and 7° that, in the absence of 
either a physical model for the y-ray emission of SFGs or a much larger 
sample of resolved galaxies, it is not possible to rule them out as a domi- 
nant contributor. 
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Figure 4: The diffuse isotropic y-ray background Black data points 
show Fermi 50-month observations ' where the horizontal bars show the 
energy bin and the vertical bars the 1 o uncertainty limit; the thick blue 
line shows our prediction for the total background due to SFGs; thin solid 
lines in other colours show the fractional contribution to this total from 7° 
decay (CR ions; green line), leptonic emission processes (CR electrons and 
positrons; orange line), and yy scattering and the resulting pair production 
cascade (red line). Broken thin blue lines show the predicted background 
if we turn off parts of our model: the blue dotted line shows the spectrum 
we would obtain if we set fcai = 1 for all galaxies at all energies, while 
the blue dashed line shows the background we would see in the absence of 


yy opacity. 


We conclude by pointing out that the methodology we have introduced 
can also be applied to predict luminosity functions and background contri- 
butions from SFGs at other wavelengths and in other messengers driven by 
CRs. Most immediately, observations by the upcoming Cherenkov Tele- 
scope Array ” and Large High Altitude Air Shower Observatory *° should 
both extend the population of y-ray-detected SFGs and push existing de- 
tections to substantially higher energies. Our model makes clear predic- 
tions for both source counts and spectral shapes that can be tested against 
these data. In the longer term, application of this model to neutrinos will 
yield predictions that will be testable by IceCube and other neutrino ob- 
servatories (see the Supplementary Information and Extended Data Figure 
3 for further information), and application to synchrotron emission from 
CR electrons can be used to make predictions for the radio sky that will be 
testable with the Square Kilometer Array (SKA) and other next-generation 
radio telescopes. Moreover, because the basis of these predictions is a co- 
herent physical model, rather than just empirical scalings, these predictions 
can all be made self-consistently. 

Supplementary Information is available for this paper. 
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Methods 


Here we describe our methods to compute y-ray emission from a single 
SFG due to both CR ions and leptons, to determine the flux received at 
Earth from that galaxy, and to apply these models to the CANDELS sam- 
ple, as well as the details of the Monte Carlo estimation for low redshift 
source counts. 


y-ray emission model for CR ions 


In our model, the total rate of y-ray emission per unit energy from a SFG 
is the sum of an ionic component and a leptonic component, dN /d Ey = 


dN.,/dE,| + dN,/dE, . We compute the ionic component as 
ion lepton 
dN, 22 1 doy dNion 
= Eion ca. Eion dEion . (1 
dEy ion " E dEy ( )| f i ( ) dEion ( ) 


Here dNion /dEion is the rate per unit energy at which supernovae injec- 
tion CR ions of energy Eion into the galaxy, opp = 40 mbarn is the mean 
proton-proton inelastic cross-section, do,/dEy(Fion) is the differential 
cross section for production of y-rays of energy E by CR ions of en- 
ergy Eion, and feail Eion) is the calorimetry fraction for CR ions of en- 
ergy Eion. We take do, /dE, (Eion) from the parameterised model of 
Ref. ?!. We compute dNion /dE;ion from the galactic star formation rate 
M. by assuming that stars form with a Chabrier initial mass function %, 
which gives the distribution of masses for newly-formed stars, and that 
stars with initial mass of 8 — 50 Mo, where Mọ is the mass of the Sun, 
end their lives as supernovae **. Each supernova injects 10°° erg of en- 
ergy in CR ions "+", distributed in energy for CR energies Fion > Mpc 
as dNion/dEion = (o M, (Pion /po) 74 dpion /dEion exp (— Eion /Ecut), 
where po = 1 GeV/c, the cutoff energy Ecut = 108 GeV, and the spec- 
tral index q = 2.2 '°?°. The exact choice of the cutoff energy above 
Eion ~ PeV makes no practical difference because the injection spectral 
index q > 2, so only a small fraction of the total CR energy is injected at 
= PeV energies regardless of Ecut, and any CRs that are injected at such 
high energies produce photons that we do not observe due to yy opac- 
ity. The normalisation factor ¢ that corresponds to our choice of initial 
mass function, supernova mass range, and CR energy per supernova is 
= 7.15 x 10% s~' GeV~* Mo" yr. 

The only remaining unknown in[Equation T]is the calorimetry fraction 
feat (Eion), which we compute from the recent model of Ref. 11 The basic 
premise of the model is that, in the neutral phase that dominates the mass 
of the ISM and thus the set of available targets for y-ray production, CR 
transport is primarily by streaming along magnetic field lines. However, 
this yields approximately diffusive transport when averaged over scales 
comparable to or larger than the coherence length of the magnetic field, 
with a diffusion coefficient D ~ Vsthg /Mi, where Vss is the CR stream- 
ing speed, hg is the gas scale height, and M4 is the Alfvén Mach number 
of the turbulence. For diffusive transport with losses in a disc geometry, the 
calorimetry fraction is given by (using the favoured parameters of Ref. !') 


9 16 A 
oF, (2.32) c] 
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where oF is the generalised hypergeometric function and Teff is the di- 
mensionless effective optical depth of the ISM, given by 


Opp Npp Ug hg C. 
2 Do Hp MH 


(3) 


Tef = 


Here npp = 0.5 is the elasticity of pp collisions, Xg is the gas surface 


density of the galactic disc, c is the speed of light, Do is the diffusion 5 


coefficient at the galactic midplane, 4p = 1.17 is the number density of 
nucleons per proton, and my = 1.67 x 10~*“ g is the mass of a hydrogen 
atom. 

To evaluate the calorimetry fraction for a CR of energy Fion, we must 
therefore determine the midplane diffusion coefficient Do for CRs of that 
energy, which in turn depends on the streaming speed Vss. This speed 
is dictated by the balance between excitation of the streaming instability 
and dissipation of the instability by ion-neutral damping, the dominant 
dissipation mechanism in the weakly-ionised neutral ISM. Balancing these 
two effects yields a CR proton streaming velocity 


ya x Ma cp?” 
4m/2 C e ura pi yatt 


where Va; is the ion Alfvén speed, ya = 4.9 x 1018 cm? gt is the ion- 
neutral drag coefficient, x is the ionised mass fraction, p = U,/2h, is the 
midplane mass density of the ISM, C is the midplane number density of 
CRs, e is the elementary charge, ui is the velocity dispersion of Alfvén 
modes in the ISM at the outer scale of the turbulence, ju; is the atomic mass 
of the dominant ion species, q is the index of the CR energy distribution, 
and y = Eion/mpc? is the Lorentz factor of the CR. Since y-ray produc- 
tion in our model is dominated by galaxies with high star formation rates 
and gas surface densities within which i) the ISM is molecule-dominated, 
ii) the main ionised species is C*, and iii) the magnetic field is set by a 
turbulent dynamo, we adopt the fiducial parameters of Ref. !! appropriate 
for such galaxies. Specifically, we take x = 10-4, i = 12, Ma = 2, 
Vai = ura / x"? Ma, and uta = Oq/V2; where gg is the gas velocity 
dispersion of the galaxy. We also adopt q = 2.2, consistent with our as- 
sumed injection spectrum. We have chosen this set of parameters without 
any fine tuning, by selecting values that are generally accepted as being the 
most appropriate for the type of source that dominates the emission. How- 
ever, we explore the parameter space in the Supplementary Information 
and show a selection of results in Extended Data Figure 1. 

At this point we have specified all the ingredients required to compute 
feai(Fion) for a galaxy of known Ny, og, and hg, save one: C, the CR 
number density. We estimate this as follows: consistent with our discus- 
sion above, for a galaxy with star formation rate per unit area ,., the CR 
ion injection rate per unit area is 


dNion : i ion “4 Eion — Ei 
= o5. f P. e Eion /Ecut dEion 
dA Mp c2 Po Pion 


Vet & min |e. Vai (1 | (4) 


(5) 
2.61 x 10° 3, s7 Myr. 
The CR number density at the midplane is then given by 
toss dNion 
Cw fee (Gin) ©) 


where tioss is the CR loss time. This is given by tioss = 1/ (tzi + tish 
where the timescale for losses in inelastic hadronic collisions is tel = 
1/ (poppNppC/Mpmu) and the diffusive escape time tait = h? Do‘. For 
the systems with high gas densities and high star formation rate that dom- 
inate y-ray production, the loss time is generally dominated by collisional 
losses. Conversely, for systems forming stars more sedately and with lower 
density environments, it is generally determined by the diffusive escape 
time. 

As a final note, we point out that the model of Ref. |! applies at 
CR energies up to tens of TeV in galaxies whose interstellar media are 
mainly molecular (most galaxies with star formation rates above a few So- 
lar masses per year), but may break down above tens of GeV in galaxies 
where the gas is mostly atomic **. Thus the model might not predict the 
correct degree of calorimetry at > 100 GeV energies in dominantly atomic 
galaxies. As shown in Extended Data Figure 2, however, low star forma- 
tion rate galaxies make only a small contribution to the background, and 
thus a possible error in them will have minimal effects on the final result. 


We illustrate the behaviour of feai (Eion) over a range of gas surface 
densities and redshifts as applied to the CANDELS sample (see below) in 
Extended Data Figures 4 and 5. 


y-ray emission model for CR leptons 


CR electrons and positrons (since both behave identically, for brevity we 
just write electrons, but everything that follows should be understood 
as applying to a mix of both) are either injected directly by diffusive 
shock acceleration in supernova remnants at the same time as CR ions 
(primary production), or appear in the decay of charged pions m™ pro- 
duced by collisions of CR ions with the ISM (secondary production). 
We assume the former carry a total energy equal to 2% of supernova 
kinetic energy '®, and have the same injection spectrum as the CR ions, 
dNe/dEe|1 x Deo 1(Ee [peje e/ Bente, with q = 2.2, but a lower spec- 
tral cut-off energy at Ecut,. = 100 TeV 35. For the latter we follow 
Refs. °°’: we first compute the rate at which CR ions produce pions of 
energy Er, 


dN, nnc dNion 
dE, = Ko B Opp (Eion) dB, 1°! (Eion) , (7) 
where nu = Yg/(2/tpmuh,) is the ISM number density, c the speed 


of light, 8 is the CR velocity divided by c, Ex = Kx(E — mpc), and 
Kx = 0.17 is the fraction of energy transferred from the CR to the pion. 
Then the rate at which these pions produce secondary electrons is 


i : 
dN, ( Ee \ dx 
=2 ey a 
3 e Lig) aE a z?’ 
where f (2) (x) is a dimensionless fititng function given in Ref. *°. 
Vu 


Thus the total CR electron injection rate is dNe /dEe = dÑe /dEe|ı + 
dNe/dEe|2. 

The electrons are subject to four dominant loss processes: collisional 
ionisation, synchrotron radiation, bremsstrahlung, and inverse Compton 
scattering; as discussed in the main text, diffusive escape from the galaxy 
is negligible in comparison. We adopt the following parameterisations 
from the literature for the total energy loss rates dEe/dt for electrons of 
energy Ee: 


dNe 
dE. 


(8) 


dE. 9 7 2, (mec? 
= = = A l zi 
di w ZOTe ms Cc ny [ina + 3 n (S) (9) 
dE: 1 
a = Tm B? 7’ 6° (10) 
sync 
dE. o 8 2 Iny+m2-1/3 yS1 
elects M vm { ®iH(1/407)/8 yZ1 a 
dE 20 
a = arcy Uraa Y (4YEpeak/MeC”). (12) 
IC 


In these expressions, ør is the Thomson cross section, œ is the fine 
structure constant, Me is the electron mass, y = Ee / Mec is the elec- 
tron Lorentz factor, Epeak is the energy where the infrared background 
peaks (derived from the dust temperature Taust = 98 (1 + ae + 
6.9 log M. /M..** and injected as a diluted modified black body spectrum 
3 Which peaks in photon number at Epeak = 2.82 kgTaust where kp 
is the Boltzmann constant), B = Vai/V4rnaupmyax is the magnetic 
field strength, uraa is the radiation energy density (which based on empir- 
ical measurements in nearby galaxies we set equal to the magnetic energy 
density w Urad = Umag = B? /87), and 1,4 and Y are dimensionless 
numerical fitting functions; these expressions are taken from Refs. *! (ion- 
isation), a2 (synchrotron), a (bremsstrahlung), and ig (inverse Compton), 
and the definitions of the fitting functions are given in those references. 


Given these loss rates, the steady-state spectrum in the galaxy is given 
approximately by 


dN. _ dNe 
dE. = dE. toss (Ee) ; (13) 


dEe 
dt 


where thoss,i = Ee 


SI ~ 

63) 2) 
Of the four loss processes, only bremsstrahlung and inverse Comp- 
ton scattering produce Ņy-rays. We compute the resulting emission using 


expressions analogous to|Equation 1| For bremsstrahlung, 


F is the loss time due to process 2, and tioss = 


1 
is the total loss time. 


dN. y NH R dN. e 
Ey, Ee dE, 14 
dE, Lee E, a OBS ( t3 ) dE. (14) 
where ogg is the cross section for production of photons of energy E, by 
electrons of energy Ee; we take our expression for ops from Refs. 444, 
Similarly, for inverse Compton we have *°“° 
dN. ra 2  dNe r 
x BE a = d G (a, ) dEs, (15) 
dEyjio 4 Epeak JEomin Ce 7 


where G(a, T) = 2a lna + (1 + 2a) (1 — a) + (Ta)? (1 — a)/2(1 +Ta), 
T = 4EpeakEe/m2c*, and a = E,/T(Ee — E,). The total lep- 


tonic contribution to y-ray production is simply dN, /dE, = 
lepton 


dÑ, /dE, + dN,/dE, | , 
brmes IC 
We show the contribution of leptonic emission to the diffuse, isotropic 
y-ray background divided by emission mechanism, and by primary versus 


secondary, in Extended Data Figure 6. 


y-ray flux at Earth 


To obtain the total observed y-ray flux for a galaxy, we must account for 
the attenuation of y-ray photons by galactic far-infrared and extragalactic 
background light (EBL) photons. We compute the optical depth 7, due to 
the former using the model of Ref. *’, and the optical depth from the latter, 
TEBL, using the model of Ref. a83 Taking these into account, we can com- 
pute the specific photon flux dF} /dE, (i.e., the number of photons per 
unit area, time and energy) received at the Earth from a galaxy at redshift 
z as 


(1+z)? dN, 
4r d? (z) dE, 


dFy _ 
dE, 


e7 EBL (Ev; z) e777 (Ey (1+2)) 


By (142) 


(16) 


where dN,,/dE, 


is the total y-ray production from both CR ions 
E,(1+z) 


and electrons evaluated at an energy E, (1+ z), and dr (z) is the luminos- 
ity distance of the source. 

The radiation that is absorbed by the host galaxy and extragalactic pho- 
ton fields is reprocessed to lower energies in the pair-production cascade. 
We parameterise the photon spectrum from this cascade using the method 
developed in Ref. °. For the purposes of our calculation here, we include 
the effect of the cascade by adding a component to dF, /dE, with a spec- 
tral shape as computed by Ref. °, and with a normalisation such that its 
energy is equal to the integrated energy lost to photon-photon scattering. 


Application of the model to CANDELS galaxies 


We apply our model to each individual galaxy in the CANDELS sample. 
The full sample contains 34,930 galaxies, but we exclude those whose 
parameters are uncertain because they contain bright active galactic nuclei, 
have unreliable redshifts, or lack a good fit to the surface brightness profile. 


6 This leaves a sample of 22,279 galaxies. 


Our calorimetry model requires, as input, the gas surface density Ug, 
scale height hg, and velocity dispersion og, along with the surface density 
of star formation Dx. However, the CANDELS data set in Ref. 7’, that 
we use, provides only the cosmological redshift z, stellar mass M., half- 
light or effective radius Re (corrected to 5000 A according to Ref. °°), and 
total star formation rate M, for our sample galaxies. We must therefore 
estimate the gas properties from observed correlations between gas and 
stellar properties. We do so as follows. 

The half light radius Re at 5000 A serves as a first order estimate of 
how the star formation and matter are distributed throughout the galac- 
tic disc. We therefore estimate the star formation rate surface density as 
D, = M, 7 2r R? and the stellar surface density as X£, = M,./ 2nR2. We 
estimate the gas surface density from the observed correlation between 
gas, stellar, and star formation surface densities given by Ref. >: 


: —0.48 
Xg — 1910-28 Da ( dx ) (17) 


Mo pe? Mo yr! pe™? \ Mo pe-? 
Similarly, there is a strong correlation between galaxy star formation rates 
and velocity dispersions, which we use to derive og. For this purpose we 
fit the relationship using the MaNGA galaxy sample °. A powerlaw fit to 
the data obtained in this survey (Fig. 6 of Ref. 5) gives 


= 18 
km s71 (18) 
Finally, we derive the gas scale height under the assumption that the gas is 
in vertical hydrostatic equilibrium, in which case the scale height is ® 


g2 


e g 
hs = OS, FE (3) 


With these gas properties in hand, we calculate an observed spectrum 
for each galaxy in the CANDELS sample by using [Equation 16] and we 
then sum over the sample to predict the y-ray flux per unit energy per unit 
solid angle ®(£,) produced by SFGs. In practice, we compute the sum 
as: 


1 Ngbin "8, dF. , 
p (Ex) = Qs > Teorey 5 ( JE- ) 7 (20) 
j=l i=1 ij 


Here Ns = 173 arcmin? is the solid angle surveyed by CANDELS, ns, j 
is the number of surveyed galaxies in the jth redshift bin, (dF, /dE,):,; 
is the flux from the ith galaxy in this bin predicted using/Equation 16] and 
Nzbin the number of redshift bins. The factor fcorr,j is the ratio of the ex- 
pected total star formation rate in each redshift bin (based on the measured 
cosmic star formation history *4 and obtained by integrating the star for- 
mation rate density over the volume in each redshift bin) to the sum of the 
star formation rates of CANDELS galaxies in that bin; its purpose is to 
correct for the fact that, due to its limited field of view and various obser- 
vational biases, the distribution of star formation with respect to redshift in 
CANDELS does not precisely match the total star formation history of the 
Universe. We use redshift bins of size Az = 0.1, chosen to ensure that the 
number of sample galaxies in each bin is large enough that the uncertainty 
in the mean spectrum due to Poisson sampling of the galaxy population is 
small. 

For the purposes of constructing Figure 2, we also require not just the 
flux, but the total y-ray luminosity in the Fermi band. We compute this 
by integrating Ey dN, /dE (computed as the sum of Equations [1] 
and[I5p from E, = 0.1 — 100 GeV. Since this comparison also requires 
the far-infrared luminosity, we convert the star formation rate to an far- 
infrared luminosity in the 8— 1000 jum band using the relation in Refs. ’*, 
corrected to a Chabrier IMF %6; this conversion is valid for star formation 
rates > 1 Mo yr 1, which encompasses almost all of the observed sample 
to which we wish to compare. 


Monte Carlo estimation for the local population 


To estimate the observable source count distribution for Fermi LAT, as 
shown in Figure 3, we cannot use the CANDELS catalogue directly, be- 
cause CANDELS has a narrow field of view that provides very little sam- 
pling of galaxies at z < 0.1, whereas these local sources are the only ones 
Fermi LAT can resolve. We therefore use a Monte Carlo scheme to sim- 
ulate a nearby, low-redshift (z < 0.1) galaxy population that follows the 
observed distribution of star formation rates in the local Universe to ac- 
count for cosmic variance, and where the correlation between galaxy star 
formation rate and y-ray luminosity is the same as what our model predicts 
for the low redshift (z < 1.5) part of the CANDELS sample. 

The first step is to produce a sample of SFGs. To do so, we draw 
galaxies from the observed distribution of star formation rates in the local 
Universe 558, For each SFG drawn, we also draw an associated redshift in 
the range z = 0 — 0.1, with probability proportional to the co-moving vol- 
ume element. We continue drawing galaxies until the total star formation 
rate of the population we have drawn matches the integrated star formation 
rate within the volume z = 0 — 0.1 as determined from the cosmic star 
formation history **. The second step is to assign y-ray luminosities for 
these galaxies based on our model for the CANDELS galaxies. For this 
purpose, we apply our model to predict the photon luminosity integrated 
over the 1 - 100 GeV band (i.e., the number of photons per unit time emit- 
ted in this energy range) for all CANDELS galaxies with z < 1.5, and fita 
power law relationship between this luminosity and the star formation rate; 
we neglect yy opacity in this calculation, since this effect is unimportant 
for the galaxies at z < 0.1 and the energy range < 100 GeV that we are 
simulating. We then assign each of our SFGs a y-ray photon luminosity 
using this powerlaw fit, and in conjunction with the redshift, an observed 
photon flux S. 

At this point we have a sample of y-ray photon fluxes S for simulated 
z < 0.1 SFGs, which we can place in bins of S to construct a synthetic 
prediction for S?(dN/dS). We carry out 13,000 Monte Carlo trials of 
this type, and in each bin of S record the mean and the 68% and 90% 
probability intervals, which we show as the blue points and bands in Figure 
3. Our method for computing the analogous confidence intervals for the 
observations is described in the Supplementary Information. 
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Galaxy [Mpc] Mo] [kpc] [Mo yr] [pc] [km s~!] [Mo pc7?] 
Arp 220 80.95 10.63 © 3.4 0 220 2 75" 100 1 4.0" 
NGC 253 3.56 > 10.4 ©! 3.87 /0.78* 4689 7/1.73%t 130! 26 |! 2.6 |! 
M31 0.77% 10.84% 2.46 © 0.26 7° 150° 10 68 0.68 5 
NGC 4945 3.725 10.52 9 2.34 70 4.07! -* 20 72 -* 


Extended Data Table 1: Local galaxy data For each entry, we give a value followed by the reference from which that value is taken. 

x NGC 4945 is observed edge-on, so measurements of the gas scale height and gas surface density are unavailable. We derive them in the usual manner, 
as described in the Methods, using the measured gas velocity dispersion. 

+ The gas data come exclusively from the nuclear starburst region, so we give two effective radii and SFR estimates: the first is for the entire galaxy, 
and the second is for the circumnuclear disk / nuclear starburst region only. We use the former for our stellar data spectrum prediction and the latter for 
our gas prediction. 
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Extended Data Figure 1: The effect of varying model parameters The plots presented here show the result of our calculations when varying the model 
parameters as discussed in the Supplementary Information. Our fiducial choice is plotted as a solid blue line, with the dashed and dash-dotted lines 
showing the spectrum for the upper and lower limits respectively of the varied parameter. The black points correspond to the Fermi data as in Figure 
4. Plot a shows Ma plotted for reasonable values of 1.6 and 2.3, and extremal values of 1.1 and 3.0; b the ionisation fraction x for values of 107? and 
1076; c the injection index q for values 2.1 and 2.3; and finally d the conversion fraction of supernova energy to CR electrons for values of 1% and 
3%, which is equivalent to 10% and 30% of the total energy injected in all cosmic ray species. Note that varying the total CR energy budget results in a 
trivial scaling of the result by the same fraction, and thus is not shown. 
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Extended Data Figure 2: The contribution of SFGs in the M, -z plane The contribution of SFGs to the total y-ray spectrum at selected energies in 
the star formation rate (M. x), redshift (z) plane. Coloured pixels show the fractional contribution (as indicated in the colourbar) from galaxies in each 
bin of M+ and z to the diffuse isotropic y-ray background at the indicated energy; a fractional contribution of unity corresponds to that pixel producing 
all of the background, with no contribution from galaxies outside the pixel. Grey points show individual CANDELS galaxies in regions of M., and z 
that contribute < 107° of the total. Flanking histograms show the fractional contribution binned in one dimension — M.: (right) and z (top). We see 
that the background at lower energies is dominated by emission from galaxies on the high side of the star forming main sequence at z ~ 1 — 2, while 
at high energies it is dominated by the brightest systems at low redshift. 
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Extended Data Figure 3: The diffuse isotropic y-ray and neutrino backgrounds The blue line and black points show the model-predicted and 
observed y-ray background, and are identical to those shown in Figure 4. The red lines show our model prediction for the neutrino background (single 


flavour) with Ecut = 100 PeV (solid line) and Ecut = 1 PeV (dashed line), computed as described in the Supplementary Information. We assume 


a neutrino flavour ratio at the detector of (ve : Vu : Vr) = (1: 1:1). The red filled band shows a power law fit ”* to the single flavour astrophysical 
neutrino background with the 90% likelihood limit, as measured by IceCube, which is also shown as grey points, where the horizontal bars show the 


energy bin and the vertical bars the 1 o uncertainty limit. 
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Extended Data Figure 4: Cosmic ray calorimetry in the E - X plane Mean calorimetry fraction feai( Æ) in the surface gas density ©, cosmic ray 
energy E plane, binned in redshift intervals. This figure is constructed by deriving the gas surface density and energy dependent calorimetry fraction 
for each galaxy in the CANDELS sample using our model. The colour of each pixel gives the mean calorimetry fraction of all the galaxies within 
that particular range of 4g, E, and redshift. The horizontal white stripes correspond to ranges of Xg into which no CANDELS galaxies fall for the 
corresponding redshift range. Several physical processes contribute to the behaviour visible in the plot. At low X4, galaxies have low feal at all energies 
E because there are few targets for hadronic collisions with CRs. As %g increases, the increased ISM density results in efficient calorimetry and 
conversion of CR energy into y-rays for low CR energies; however, at higher energies the CR number density is low, yielding a high CR streaming 
velocity and rapid escape, resulting in low fcai. As Xg increases further, the increasing density results in the streaming instability being suppressed 
efficiently by ion-neutral damping towards lower energies, reducing the calorimetry fraction further. Finally, at the highest XJ, the streaming instability 
is suppressed completely by ion-neutral damping, but streaming is still limited to the speed of light. Consequently, increasing i, further only results in 
increased collisions, and thus a higher calorimetry fraction. 
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Extended Data Figure 5: Cosmic ray calorimetry in the z - X plane Mean calorimetry fraction in the surface gas density (%4), redshift (z) plane 
at CR energies Æ = 1 GeV, 10 GeV, 1 TeV and 10 TeV. To construct this figure, for each CANDELS sample galaxy, we apply our model to compute 
Ng and fcai( E) at the indicated energies. The colour indicates the average fcai() value computed over bins of (z, 4), while contours indicate the 
density of the CANDELS sample in this plane. Note that the non-monotonic behaviour of fea( Æ) with £g that is most prominently visible in the 1 
TeV panel is expected, for the reasons explained in the caption of Extended Data Figure 4. 
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Extended Data Figure 6: Contributions to the diffuse isotropic y-ray background The blue line and black points show the model-predicted and 
observed y-ray background, and are identical to those shown in Figure 4. The green line shows the contribution from 7° decay, the olive lines the 
contribution from bremsstrahlung emission, and the cyan lines the contribution from the inverse Compton emission. In both cases, dashed lines show 


the spectrum produced by primary CR electrons and the dash-dotted lines the spectrum from secondary electrons and positrons. The red line shows the 
contributions from the EBL cascade. 
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Supplementary Information 


1 Confidence intervals for source count distributions 


Calculation of confidence intervals on S?(dN/dS) (as shown in Figure 3) for the observed sources (both the Fermi-observed SFGs and our model- 
predicted CANDELS SFGs) is non-trivial, because both surveys cover a fraction f < 1 of the sky, and the number of sources per bin for at least some 
bins of S' is small, so we cannot compute the uncertainty by assuming that we are in the large N limit. To perform the calculation, we assume that 
the SFG population follows a Poisson distribution on the sky (i.e., we are in the cosmological isotropic limit), so if the entire sky contains Niot SFGs 
within some bin of photon flux, the probability that Nops will be found within the observable region can be written 


(f Neo) Vor e7 f Ntot 
Nors! 


P (Nobs|Ntot) = (21) 


We wish to solve the inverse problem, i.e., given an observed number Nops, what is the probability distribution of Ntot? The answer is given by 
Bayes’s Theorem, which requires 


P (Ntot|Nobs) = P (Nobs|Netot (22) 


where P (Ntot|Nops) is the posterior probability, P (Ntot) is the prior probability, and P (Vops) is a normalisation factor. We adopt a flat prior 
P (Mtot) « 1, so we can then write 


P (Niot| Nove) = MN Nope TNs, (23) 


where M is a normalisation constant. For e7% < 1, which is always the case since 0 < f < 1, the value of M required to guarantee that 
D No P (Neot | Nobs) = lis 


N = — r (24) 


where Lis (z) is the polylogarithm of order s. 


To compute the confidence interval we require the cumulative distribution function. In the discrete case this is given by calculating the probability 
that Ntot < N, which is 


N-1 
PORN SSN See (25) 
n=0 
= 1-vS> nNovse—fn (26) 
n=N 
= 1 -N5 (i + N)Novs e7S CHN) (27) 
i=0 


et & (ef, — Nobs, N) 
Li- No»s (ef) 


1 (28) 


where ® (z, s, a) is the Lerch Phi function (sometimes also referred to as the Lerch Zeta function). To obtain a particular percentile p in the range 0 to 
1, we simply use the continuous forms of the polylogarithm and the Lerch Phi functions, set p = P (Ntot < N) and invert the problem numerically to 
find the appropriate value for N; for the purposes of Figure 3, we are interested in the 90% confidence interval, so we take p = 0.05 and p = 0.95. For 
the special case Nops = 0, the result simplifies to 


een 1 
p=l-e (z (e-F) + 1) ; (29) 


which we can invert numerically for p = 0.9 to obtain the 90% confidence upper limit. 

In order to use the result we have just derived, we require a value for f. For the CANDELS data points, this is straightforward: our data come from 
the GOODS-S field, which has an area of 173 arcmin”, corresponding to f = 1.16 x 107°. Assigning a value of f to the Fermi data is more complex: 
Fermi LAT surveys the entire sky, but it cannot detect faint sources, such as SFGs, that are too close to the Galactic plane because they are hidden by the 
Galactic diffuse foreground. As a result, the effective survey area depends at least somewhat on the flux and spectral shape of the target SFG — brighter 
and harder sources can be detected closer to the plane than fainter and softer ones. Capturing this effect in detail would require extensive testing of the 
Fermi reduction pipeline using artificial sources, which is beyond the scope of this work. For the purposes of computing the confidence intervals shown 
in Figure 3, we ignore this complexity, and roughly estimate that SFGs are undetectable within 15° degrees of the Galactic plane, which corresponds to 


approximately f = 0.7. 15 


2 Sensitivity of the result to model parameters 


Here we investigate the sensitivity of our results to our choice of model parameters. The set of parameters we have adopted to derive the key result of 
this study were deliberately chosen to be approximately the consensus value. We purposefully chose not to fine-tune our parameters to obtain a best 
fit, nor, more importantly, have we had to adopt extremal values of any parameter to force the calculation in the direction of making the contribution 
of SFGs dominant. However, it is still of interest to explore the parameter space within a reasonable range of variation. The key tunable parameters in 
our model are: (1) the total energy per unit mass of stars formed that is ultimately injected in CRs, (2) the Alfvén Mach number Ma, (3) the ionisation 
fraction x, (4) the CR injection spectral index q, and (5) the ratio of CR energy injected into primary leptons to that injected into primary protons. We 
discuss each of these in turn. 

The total CR energy budget is constrained by observations of both individual SN remnants '° and by observation of the total y-ray luminosity of 
local starburst galaxies, which are generally thought to be fully calorimetric '°”>**. These observations constrain this energy budget to be within a 
factor of 2 of our fiducial choice. Varying this parameter within that range would increase or decrease our predicted background from SFGs by the 
same factor, while leaving the overall shape of the spectrum unchanged. 

The Alfvén Mach number Ma cannot be measured directly in external galaxies, but is determined by the ratio of the magnetic ug and kinetic Ukin 
densities in the flow, Ma = (3up/ Ukin) °° 11 Which in turn is well constrained by dynamo theory and simulations 117475 Galactic magnetic fields are 
driven by both turbulent and aQ dynamos. The growth timescale for the latter is the orbital period, and for the former is the eddy turnover time, which 
for a galaxy that is marginally stable against gravity is also comparable to the orbital period ”°; since essentially all galaxies are many orbital periods old, 
we can assume that galactic dynamos have reached saturation. For the turbulent dynamo driven by supersonic motion, the saturation level of ug /Ukin 
is determined by the Reynolds and magnetic Prandtl numbers. Both of these dimensionless numbers are large in star-forming galaxies (Re = 10° 76 
and Pr ~ hg/lap = 100, where hg is the gas scale height and Jp is the length scale at which magnetic fields decouple from the gas due to ambipolar 
diffusion Ny. In this regime, Ref. 75 find up /Uxin = 0.040 — 0.064, which corresponds to Alfvén Mach numbers of 2.9 and 2.3. This provides an 
absolute upper limit on M4, expected from the turbulent dynamo alone. The aQ dynamo will further enhance the field strength by creating a coherent 
component on top of the turbulent one. Simulations in Ref. ™ find that the aQ dynamo in galactic discs saturates at Boon/Bturb = 0.27 — 0.42, 
corresponding to a factor of 1.27? — 1.42? increase in ug compared to our estimate for the turbulent dynamo only. If we take the minimum value for the 
turbulent energy from Ref. ” (0.04) and the minimum af amplification factor from Ref. ”* (1.277), this gives Ma = 2.27, while the maximal values 
0.064 and 1.42? give Ma = 1.61. Our fiducial choice of Ma = 2 is simply the middle of this fairly narrow range of reasonable values. We show the 
results of varying Ma in Panel a of Extended Data Figure 1, where we find that lowering Ma to 1.6 yields a slightly fainter background, while raising 
it to Ma = 2.3 makes the predicted background somewhat brighter; the overall shape is largely unchanged in either case. For illustrative purposes we 
also include the extremal values of 1.1 and 3.0, corresponding to near-equipartition and very sub-equipartition field strengths. 

The ionisation fraction x has a larger plausible range: in galaxies with predominantly atomic interstellar media such as the Milky Way, it reaches 
x =~ 1077, while in extreme starbursts with very high densities it might reach as low as 10~°, indicating a 2 dex range around our fiducial choice ''”. 
However, this parameter also enters the problem to the 1/4 power. We show the effects of varying x in Panel b of Extended Data Figure 1; it is clear 
from this figure that, despite its larger range of variation, varying x within its plausible range has a smaller effect than varying Ma. 

The injection spectral index q is constrained by observations of local SN remnants, which require that it lie in the range q ~ 2.1 — 2.3 1°. 
We explore this range of variation in Panel ¢ of Extended Data Figure 1, which shows that changing the injected spectral index induces, as might be 
expected, a similar change in the spectral slope of the predicted background. However, given the uncertainties, it is not clear exactly which spectral 
index in this range would give the best agreement with the data, and we note that all of these models lie far from what would be expected for full 
calorimetry, as is clear from comparing Extended Data Figure 1 to the line for full calorimetry in Figure 4 of the main text; even for q = 2.1 or 2.3, the 
overall spectral slope is determined mostly by variation of fcai(/), not by the choice of injection spectrum. 

Finally, our fiducial choice of ratio of primary electrons to protons is motivated by what is required to reproduce the FIR-radio correlation !°. 
However, given various uncertainties, this could plausibly change at the factor of ~ 2 level. We show the effects of varying the fraction of SN energy 
that goes into primary CR electrons (while holding CR ions constant) in Panel d of Extended Data Figure 1. As expected, the effects are minimal except 
below = 1 GeV in energy, and even there are small, since reducing the energy in primary electrons of course does not affect emission from secondaries, 
which are equally important. 

The overall message of Extended Data Figure 1 is that the conclusion that SFGs contribute at least ~ 50% of the total diffuse y-ray background 
seems inescapable, even if we choose extreme values for model parameters. However, there is also another, somewhat stronger conclusion to draw: 
within the reasonable parameter space, the spectral shape we derive, which is dictated primarily by our model for feaı( Æ), matches the spectral shape 
of data well over a > 3 decade range in y-ray energy. This means that, if we were to adopt extreme values of parameters such that SFGs contribute only 
= 50% of the background, leaving the remainder to be made up by some other unknown source population, matching the observed spectrum would 
require implausible fine-tuning: this other source population would have to produce a spectral shape and magnitude nearly identical to that of SFGs, 
over the entire energy range from ~ 0.3 GeV to ~ 1 TeV. 

The one place where there is some minor tension between the spectral shape predicted by our model and that observed is near the EBL-induced 
cutoff of the spectrum between ~100 GeV and ~1 TeV, where our model falls slightly below the data (though the data themselves are uncertain in this 
energy range due to the difficulties of background subtraction). One possible explanation is that the model prediction in this energy range is sensitive 
to our chosen functional form for the EBL optical depth, and could conceivably be fit better by alternative models. Moreover, due to the EBL, the 
background at these energies is dominated by low-redshift, hard y-ray sources (see Extended Data Figure 2), which are poorly sampled by the deep, 
narrow field of view in CANDELS. We have partly accounted for this by correcting the CANDELS star formation density using the measured star 
formation history, but a better solution, which we leave for future work, would be to supplement CANDELS by a wide-field, low-redshift galaxy 
survey. 


3 Comparison to earlier work 


Our conclusion that emission from SFGs dominates the diffuse, isotropic y-ray differs from some earlier work. It is therefore important to examine 
the precise reasons why this is the case. One contributing factor is that eeplier models were forced to adopt single power laws for the emitted y-ray 


spectrum in different classes of galaxies “*”. By contrast, Figure 1 demonstrates that none of the four nearby resolved galaxies shown have spectra 
that are well described by a y-ray spectrum in the form of a pure power law over the energy range from Æ, = 1 — 1000 GeV; our model correctly 
captures this behaviour, but earlier pure power law models did not. Similarly, we calculate fecal as a function of energy directly, rather than relying on an 
empirical FIR-y correlation, and our calorimetry fractions are on average larger than those implicitly assumed in earlier works. This is because many 
of the lower estimates for the contribution from SFGs to the y-ray background rely on a FIR-y correlation derived from early Fermi detections of < 10 
individually-resolved SFGs ’ that yields somewhat lower y-ray luminosities than more recent fits using a larger (but still small) sample of SFGs ”, and 
with which our model agrees (Figure 2). Thus the reason we find that SFGs can produce the full background, whereas earlier models could not, is that 
our model predicts y-ray emission that is both somewhat brighter and has a more complex spectral shape than the values adopted in earlier work. 

Likewise, earlier claims that a variety of other source classes dominate the diffuse, isotropic background have also relied on extrapolated empirical 
correlations with large uncertainties. For instance Ref. ® estimates the contribution from misaligned active galactic nuclei using a radio-y correlation 
derived from a sample of 16 resolved objects, coupled with a radio luminosity function extrapolated to redshifts considerably higher than those well- 
sampled by observations *'. By contrast, our assignment of y-ray luminosities to SFGs is based on a physical model that agrees with local observations, 
and the CANDELS catalogue from which we draw our distribution of SFG properties has very good completeness over the range of redshift and star 
formation rate that dominates production of the diffuse, isotropic y-ray background (see Extended Data Figure 2). 


4 Neutrinos 


In addition to electrons and positrons, the decay of 7~ also produces neutrinos, which are of particular interest as they propagate largely unhindered 
from the source to the observer. Our goal here is to compute the all-species neutrino flux due to SFGs, so that we may compare to the astrophysical 
neutrino background measured by IceCube ”. 

The relationship between the y-ray and neutrino spectra is approximately given by E;F, (E, = E,/2) = (3/2)E>F,(E,) *°. However, we 
compute the neutrino flux from the charged pion decay in our sample galaxies using the more detailed method in Ref. ial | 

Charged pion decay produces neutrinos in two steps: the initial decay of the pion creates a muon and a muon neutrino, and then the muon decays, 
yielding an electron, an electron neutrino, and a second muon neutrino (where we do not distinguish between particles and anti-particles). The all-flavour 
neutrino spectrum is then a sum over the energy distributions of all three neutrinos produced in this chain, given by 


dN, 1 dN, (Ev,\ dx 2 f* dN, (E \ dz 
m =f (I. (+ ho (°)) dEr ( x ) oo J dEr ( T ) r’ (30) 


where A = 1 — (m,/ mr)’, x = E,/E;,, the second integral accounts for the muon neutrinos produced in the initial charged pion decay, and the 
functions fy, and fi) describe the energy distributions for the electron and muon neutrinos produced by decay of the secondary muon, respectively; 
ù 


we take them from Equations 40 and 36 of Ref. °°. The ratio of neutrino flavours at the source is (Ve : Vu : Vr) = (1 : 2: 0). However, neutrino 
oscillations will bring this to an even (ve : Vu : Vr) = (1: 1: 1) for an observer at Earth. 

[Equation 30ļis the analogue to Equation 1 of the main text for y-rays, and we can compute the resulting specific neutrino flux for each galaxy from 
Equation 2 of the main text simply by replacing dN, /dE., with dN, /dE, and setting the opacities tyy = TeBL = 0. We use this to calculate a 
predicted neutrino flux from each CANDELS galaxy, and we sum to compute the neutrino background due to SFGs using Equation 3 of the main text, 
exactly as we do for the y-ray background. We plot the resulting predicted neutrino spectrum in Extended Data Figure 3. 

We find that our model predicts that SFGs produce a neutrino flux that is ~ 15% of the astrophysical neutrino background, as measured by IceCube 
7 for a CR spectral cutoff energy of Ecut = 100 PeV. However, while the choice of cut has no significant effect on the y-ray spectrum (as explained in 
the main text), it does matter for the neutrino spectrum due to the high energies of the astrophysical neutrinos observable by IceCube. Consequently, we 
find that SFGs produce < 15% of the observed neutrino background if we adopt a smaller value of Ecut 7°. To illustrate this, in Extended Data Figure 3 
we show two calculations: one with our fiducial Ecut = 100 PeV, and one with a smaller Ecut = 1 PeV. The cutoffs in the neutrino spectrum shown 
in Extended Data Figure 3 are a direct result of the adopted value of Ecut. We also note that the normalisation of our predicted neutrino spectrum is 
sensitive to bright, hard neutrino sources at low redshift, which dominate at high energy but are poorly sampled by the small CANDELS field of view. 
This suggests that it would be worthwhile in the future to repeat this analysis using a survey of SFGs that is wider but shallower than CANDELS. A 
further consideration relates to the exact form of the cosmic ray calorimetry at ultra high energies. Here we would expect cosmic rays with sufficiently 
large gyro radii to interact with and scatter off the external turbulence above the damping scale, significantly decreasing the diffusion coefficients and 
increasing calorimetry in turn. We leave an exploration of this mechanism and the injection cutoff energy for future work. 
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